/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package methods;

import core.Constants;

/**
 *
 * @author Alexey Chirkov
 */
public class KrankaNikolsona implements Method {

    private static KrankaNikolsona obj = null;
    private int K = 4;
    private int I = 4;
    private double ht = 45;
    private double hr = 1;
    private double[][] value;

    public KrankaNikolsona() {
        solved();
    }

    public static Method getInstance() {
        if (obj == null) {
            obj = new KrankaNikolsona();
        }
        return obj;
    }

    private void solved() {
        value = new double[K + 1][I + 1];
        for (int i = 0; i <= I; i++) {
            value[0][i] = Constants.getU0();
        }
        double[][] matA = new double[I + 1][I + 1];
        double[] matB = new double[I + 1];

        for (int k = 1; k <= K; k++) {
            matA[0][0] = getY1() + 4 * getY0();
            matA[0][1] = -4 * getY0() + getY3();
            matB[0] = value[k-1][1]*4*getY0()+value[k-1][0]*(getY1()-4*getY0()-getY3())+getI0(0);
            for (int i = 1; i < I; i++) {

                matA[i][i - 1] = -getY0() + (getY0() / i);
                matA[i][i] = getY1() + 2 * getY0() - (getY0() / i) + getY3();
                matA[i][i + 1] = -getY0();
                matB[i] = value[k-1][i+1]*getY0()+value[k-1][i]*(getY1()-2*getY0()+(getY0()/i)-getY3())+value[k-1][i-1]*(getY0()-(getY0()/i))+getI0(i);
            }
            matA[I][I - 1] = -2 * getY0() + getY2();
            matA[I][I] = getY1()+2*getY0()-getY2()+getY3();
            matB[I] = value[k-1][I]*(getY1()-2*getY0()+getY2()-getY3())+value[k][I-1]*(2*getY0()-getY2())+getI0(I);
            double[] matX = propaceMothod(matA, matB);
            for (int i = 0; i <= I; i++) {
                value[k][i] = matX[i];
            }
        }
    }

    private double getY0() {
        double temp = 0.5 * Constants.getK() / Constants.getC();
        return ht * temp;
    }

    private double getY1() {
        return hr * hr;
    }

    private double getY2() {
        double temp1 = Constants.getK() * ht * hr;
        double temp2 = 2 * Constants.getC() * Constants.getRadius();
        return temp1 / temp2;
    }

    private double getY3() {
        double temp1 = Constants.getK() * Constants.getAlpha();
        double temp2 = Constants.getC() * Constants.getL();
        double temp = temp1 / temp2;
        return temp * ht * hr * hr;
    }

    private double getI0(int i) {
        double temp1 = 2 * Constants.getK() * Constants.getAlpha();
        double temp2 = Constants.getC() * Constants.getL();
        double temp = temp1 / temp2;
        double sum1 = temp * Constants.getU0() * ht * hr * hr;
        temp = Constants.getBetta() / Constants.getC();
        double sum2 = temp * Constants.getIntensity(getRi(i)) * ht * hr * hr;
        return sum1 + sum2;
    }

    @Override
    public double getValueAt(int k, int i) {
        return value[k][i];
    }

    @Override
    public void setK(int k) {
        this.K = k;
        this.ht = Constants.getTimeEnd() / k;
        solved();
    }

    @Override
    public int getK() {
        return K;
    }

    @Override
    public void setI(int i) {
        this.I = i;
        this.hr = Constants.getRadius() / i;
        solved();
    }

    @Override
    public int getI() {
        return I;
    }

    @Override
    public double getRi(int i) {
        return i * hr;
    }

    @Override
    public double getTk(int k) {
        return k * ht;
    }

    private double[] propaceMothod(double a[][], double[] b) {
        int N = a.length;
        double m;
        for (int i = 1; i < N; i++) {
            m = a[i][i - 1] / a[i - 1][i - 1];

            a[i][i - 1] = 0;
            a[i][i] = a[i][i] - m * a[i - 1][i];
            if (i != N - 1) {
                a[i][i + 1] = a[i][i + 1] - m * a[i - 1][i + 1];
            }
            b[i] = b[i] - m * b[i - 1];
        }
        double[] x = new double[N];
        for (int i = N - 1; i >= 0; i--) {
            if (i == N - 1) {
                x[i] = b[i] / a[i][i];
            } else {
                x[i] = (b[i] - a[i][i + 1] * x[i + 1]) / a[i][i];
            }
        }
        return x;
    }
}
